# 大型钢框架结构 - 10层6×6跨（约360节点，1000+自由度）
# 基于原始模型扩展，适用于GPU求解器测试

wipe
model basic -ndm 3 -ndf 6

# 材料定义（国际单位制：m-N-kg）
# 混凝土材料 Concrete02
uniaxialMaterial Concrete02 1 -26.8e6 -0.002 -10.0e6 -0.006 0.1 2.2e6 44.0e9

# 钢材材料 Steel01 (Q345钢)
uniaxialMaterial Steel01 2 345.0e6 206.0e9 0.01

# 几何参数（国际单位制：米）
set L1 4.0;        # X向跨度 4m
set L2 3.5;        # Y向跨度 3.5m  
set H  3.0;        # 层高 3m
set numFloorsX 6;  # X向跨数（7个节点）
set numFloorsY 6;  # Y向跨数（7个节点）
set numStories 100; # 总层数

# 输出模型信息
set totalNodes [expr ($numFloorsX + 1) * ($numFloorsY + 1) * ($numStories + 1)]
set totalDOF [expr $totalNodes * 6]
puts "模型规模："
puts "楼层数: $numStories"
puts "X向跨数: $numFloorsX (节点数: [expr $numFloorsX + 1])"
puts "Y向跨数: $numFloorsY (节点数: [expr $numFloorsY + 1])"
puts "总节点数: $totalNodes"
puts "总自由度: $totalDOF"

# 节点生成函数
proc generateNodes {numX numY numZ L1 L2 H} {
    set nodeTag 1
    for {set k 0} {$k <= $numZ} {incr k} {
        for {set j 0} {$j <= $numY} {incr j} {
            for {set i 0} {$i <= $numX} {incr i} {
                set x [expr $i * $L1]
                set y [expr $j * $L2]
                set z [expr $k * $H]
                node $nodeTag $x $y $z
                incr nodeTag
            }
        }
    }
}

# 节点编号计算函数（i,j,k对应x,y,z方向索引）
proc getNodeTag {i j k numX numY} {
    return [expr $k * ($numX + 1) * ($numY + 1) + $j * ($numX + 1) + $i + 1]
}

# 生成所有节点
puts "正在生成节点..."
generateNodes $numFloorsX $numFloorsY $numStories $L1 $L2 $H

# 底层全固定约束
puts "设置边界约束..."
set nodesPerFloor [expr ($numFloorsX + 1) * ($numFloorsY + 1)]
for {set i 1} {$i <= $nodesPerFloor} {incr i} {
    fix $i 1 1 1 1 1 1
}

# 几何变换定义
geomTransf Linear 1 1 0 0;  # 柱单元变换
geomTransf Linear 2 0 1 0;  # X方向梁单元变换
geomTransf Linear 3 0 0 1;  # Y方向梁单元变换

# 截面参数（保持原模型参数）
# 梁截面 H200×600×20×20
set bf_beam 0.200
set h_beam 0.600
set tf_beam 0.020
set tw_beam 0.020
set A_beam [expr 2*$bf_beam*$tf_beam + ($h_beam-2*$tf_beam)*$tw_beam]
set Iz_beam [expr ($bf_beam*pow($h_beam,3) - ($bf_beam-$tw_beam)*pow($h_beam-2*$tf_beam,3))/12]
set Iy_beam [expr (2*$tf_beam*pow($bf_beam,3) + ($h_beam-2*$tf_beam)*pow($tw_beam,3))/12]
set J_beam [expr 2*$bf_beam*pow($tf_beam,3)/3 + ($h_beam-2*$tf_beam)*pow($tw_beam,3)/3]

# 柱截面 H600×600×20×20
set bf_col 0.600
set h_col 0.600
set tf_col 0.020
set tw_col 0.020
set A_col [expr 2*$bf_col*$tf_col + ($h_col-2*$tf_col)*$tw_col]
set Iz_col [expr ($bf_col*pow($h_col,3) - ($bf_col-$tw_col)*pow($h_col-2*$tf_col,3))/12]
set Iy_col [expr (2*$tf_col*pow($bf_col,3) + ($h_col-2*$tf_col)*pow($tw_col,3))/12]
set J_col [expr 2*$bf_col*pow($tf_col,3)/3 + ($h_col-2*$tf_col)*pow($tw_col,3)/3]

# 材料常数
set E_steel 206.0e9
set G_steel 79.0e9

puts "截面特性："
puts "梁截面: A=[format %.6f $A_beam]m², Iz=[format %.3e $Iz_beam]m⁴"
puts "柱截面: A=[format %.6f $A_col]m², Iz=[format %.3e $Iz_col]m⁴"

# 生成柱单元
puts "正在生成柱单元..."
set elemTag 10000
for {set k 0} {$k < $numStories} {incr k} {
    for {set j 0} {$j <= $numFloorsY} {incr j} {
        for {set i 0} {$i <= $numFloorsX} {incr i} {
            set node1 [getNodeTag $i $j $k $numFloorsX $numFloorsY]
            set node2 [getNodeTag $i $j [expr $k+1] $numFloorsX $numFloorsY]
            element elasticBeamColumn $elemTag $node1 $node2 $A_col $E_steel $G_steel $J_col $Iy_col $Iz_col 1
            incr elemTag
        }
    }
}

# 生成X向梁单元
puts "正在生成X向梁单元..."
for {set k 1} {$k <= $numStories} {incr k} {
    for {set j 0} {$j <= $numFloorsY} {incr j} {
        for {set i 0} {$i < $numFloorsX} {incr i} {
            set node1 [getNodeTag $i $j $k $numFloorsX $numFloorsY]
            set node2 [getNodeTag [expr $i+1] $j $k $numFloorsX $numFloorsY]
            element elasticBeamColumn $elemTag $node1 $node2 $A_beam $E_steel $G_steel $J_beam $Iy_beam $Iz_beam 2
            incr elemTag
        }
    }
}

# 生成Y向梁单元
puts "正在生成Y向梁单元..."
for {set k 1} {$k <= $numStories} {incr k} {
    for {set j 0} {$j < $numFloorsY} {incr j} {
        for {set i 0} {$i <= $numFloorsX} {incr i} {
            set node1 [getNodeTag $i $j $k $numFloorsX $numFloorsY]
            set node2 [getNodeTag $i [expr $j+1] $k $numFloorsX $numFloorsY]
            element elasticBeamColumn $elemTag $node1 $node2 $A_beam $E_steel $G_steel $J_beam $Iy_beam $Iz_beam 3
            incr elemTag
        }
    }
}

# 质量定义（楼板荷载分布到节点）
puts "设置质量..."
# 每个节点影响面积：L1 × L2 / 4 = 4.0 × 3.5 / 4 = 3.5m²
# 楼板荷载：6.5kN/m² = 6500N/m²
# 节点重量：6500N/m² × 3.5m² = 22750N
# 节点质量：22750N / 9.8m/s² = 2321.4kg
set nodeMass 2321.4

# 给除底层外的所有节点分配质量
for {set k 1} {$k <= $numStories} {incr k} {
    for {set j 0} {$j <= $numFloorsY} {incr j} {
        for {set i 0} {$i <= $numFloorsX} {incr i} {
            set nodeTag [getNodeTag $i $j $k $numFloorsX $numFloorsY]
            mass $nodeTag $nodeMass $nodeMass $nodeMass 0 0 0
        }
    }
}

# 静力荷载模式
set nodeLoad 22750.0
pattern Plain 1 Linear {
    for {set k 1} {$k <= $numStories} {incr k} {
        for {set j 0} {$j <= $numFloorsY} {incr j} {
            for {set i 0} {$i <= $numFloorsX} {incr i} {
                set nodeTag [getNodeTag $i $j $k $numFloorsX $numFloorsY]
                load $nodeTag 0 0 [expr -$nodeLoad] 0 0 0
            }
        }
    }
}

# 静力分析设置
constraints Transformation
numberer RCM
system BandGeneral
test NormDispIncr 1e-6 100 0 1
algorithm Newton
integrator LoadControl 0.1
analysis Static

puts "开始静力分析..."
set staticOK [analyze 10]

if {$staticOK == 0} {
    puts "静力分析成功完成"
} else {
    puts "静力分析失败"
    exit
}
loadConst -time 0.0

# 输出关键节点位移
set midX [expr $numFloorsX / 2]
set midY [expr $numFloorsY / 2]
puts "\n静力分析结果（中心线节点位移）："
for {set k 1} {$k <= $numStories} {incr k} {
    set nodeTag [getNodeTag $midX $midY $k $numFloorsX $numFloorsY]
    puts "第${k}层中心节点${nodeTag}: Z位移 = [format %.6f [nodeDisp $nodeTag 3]] m"
}

# 地震动分析准备
wipeAnalysis
remove loadPattern 1

# 地震动时程设置
set dt 0.01
timeSeries Path 1 -dt $dt -filePath "earthquake.txt" -factor 3.0
pattern UniformExcitation 2 1 -accel 1

# 瑞利阻尼（模态阻尼比5%）
set eigenValues [eigen 3]
set lambda1 [lindex $eigenValues 0]
set lambda2 [lindex $eigenValues 1]
if {$lambda1 > 0 && $lambda2 > 0} {
    set omega1 [expr sqrt($lambda1)]
    set omega2 [expr sqrt($lambda2)]
    set xi 0.05
    set a0 [expr $xi*2.0*$omega1*$omega2/($omega1+$omega2)]
    set a1 [expr $xi*2.0/($omega1+$omega2)]
    rayleigh $a0 0.0 0.0 $a1
}

# 记录关键节点响应
set topNodeTag [getNodeTag $midX $midY $numStories $numFloorsX $numFloorsY]
recorder Node -file "displacement_x_large.out" -time -node $topNodeTag -dof 1 disp
recorder Node -file "acceleration_x_large.out" -time -node $topNodeTag -dof 1 accel

# 动力分析设置
constraints Transformation
numberer RCM
# system BandGeneral
system GPUSolver
test NormDispIncr 1e-6 25 0
algorithm ExpressNewton 2 1 1.0
integrator Newmark 0.55 0.275625
analysis Transient

# 执行地震分析
set analysisTime 10.0; 
set numSteps [expr int($analysisTime/$dt)]
puts "开始地震动时程分析，共$numSteps步..."

set dynamicStartTime [clock milliseconds]
set dynamicOK [analyze $numSteps $dt]
set dynamicEndTime [clock milliseconds]
set dynamicTime [expr $dynamicEndTime - $dynamicStartTime]

if {$dynamicOK == 0} {
    puts "地震动时程分析成功完成"
} else {
    puts "地震动时程分析完成（可能有收敛问题）"
}
puts "动力分析耗时: ${dynamicTime} 毫秒 ([format %.3f [expr $dynamicTime/1000.0]] 秒)"

# # 输出最终结果
# puts "\n地震动分析最终位移结果："
# for {set k 1} {$k <= [expr min($numStories, 5)]} {incr k} {
#     set nodeTag [getNodeTag $midX $midY $k $numFloorsX $numFloorsY]
#     puts "第${k}层中心节点${nodeTag}: X位移 = [format %.6f [nodeDisp $nodeTag 1]] m"
# }

# # 模态分析
# set eigenValues [eigen 5]
# puts "\n模态分析结果："
# for {set i 0} {$i < [llength $eigenValues]} {incr i} {
#     set eigenVal [lindex $eigenValues $i]
#     if {$eigenVal > 0} {
#         set freq [expr sqrt($eigenVal)/(2*3.14159)]
#         set period [expr 1.0/$freq]
#         puts [format "第%d阶: 频率 = %.3f Hz, 周期 = %.3f 秒" [expr $i+1] $freq $period]
#     }
# }

# puts "\n大型模型创建完成！"
# puts "总节点数: $totalNodes"
# puts "总自由度: $totalDOF"
# puts "总单元数: [expr $elemTag - 10000]"
